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We consider the problem of controlling the motion of an atom trapped in an optical cavity using 
continuous feedback. In order to realize such a scheme experimentally, one must be able to perform 
state estimation of the atomic motion in real time. While in theory this estimate may be provided 
by a stochastic master equation describing the full dynamics of the observed system, integrating 
this equation in real time is impractical. Here we derive an approximate estimation equation for 
this purpose, and use it as a drive in a feedback algorithm designed to cool the motion of the atom. 
We examine the effectiveness of such a procedure using full simulations of the cavity QED system, 
including the quantized motion of the atom in one dimension. 



PACS numbers: 02.30.Yy,32.80.Pj,42.50.-p,03.67.-a 
I. INTRODUCTION 

The theory of quantum measurement has been mired 
in controversy for the majority of its history. In a large 
part this is due to the philosophical difhculties, usually 
referred to as the "quantum measurement problem," with 
which it is associated 0|. However, in most areas of 
physics it has been generally possible to ignore quantum 
measurement theory, because it was largely irrelevant to 
real measurements that were being made by experimen- 
talists 0]- In particular, if needed at all, an ensemble pic- 
ture in which quantum measurement theory merely pre- 
dicted average statistics provided a sufficient description. 
However, experimental technology has now advanced to 
the point where repeated measurements may be made on 
a single quantum system as it evolves ^,4s«5»i£(|- As a re- 
sult, not only does quantum measurement theory become 
directly relevant to experimental physics, but we can con- 
sider testing the predictions of this theory regarding the 
effects of measurement on the single system itself. 

While the theory of quantum measurement is now 
widely accepted, some still remain unconvinced that the 
conditioned state describes what is really happening to 
an individual system when it is measured. How might we 
verify these predictions, and bring at least this part of the 
controversy to an end? One way to do this is to attempt 
to perform feedback control on an individual quantum 
system 0, ■ Such an experiment involves taking into 
account the changes in the system due to the measure- 
ment results, and using this knowledge to alter inputs to 
the system in order to control its dynamics. This would 
therefore provide a direct verification of quantum mea- 
surement theory: if the feedback algorithm works, then 
the theory is correctly predicting the behavior of the sys- 
tem; conversely, if the theory is not giving accurate pre- 



dictions, the feedback algorithm will be ineffective. 

Cavity quantum electrodynamics (CQED) is one area 
in which the motion of a single quantum system, an in- 
dividual atom in an optical cavity, can be monitored 
experimentally in real time |^ ^ IE, IE IS- Not only 
does the light emitted from the optical cavity provide 
a means to track the atomic motion, but the laser driv- 
ing the cavity provides a means to apply forces to the 
atom. This system is therefore an excellent candidate 
for implementing real-time quantum feedback control. It 
is our purpose here to examine the implementation of 
such a process in this system, and to show that by us- 
ing approximate estimation techniques, which will be es- 
sential in such experiments for the foreseeable future, 
a measurable degree of control over the atomic mo- 
tion can be realized. In particular, we will simulate a 
control algorithm designed to cool the atomic motion. 
^^is^i^l^e^^ferred^C) ^j'^j^^ "^ ^ l^OlT ll^l'^^ l'^si 

sive cooling schemes which have also been the focus of 
recent theoretical UllllllilllllllliliSlillliland 
experimental |43j | work. 

In general, feedback control of the atomic motional 
state requires the real-time estimation of the quantum 
state from the continuous measurement, thus provid- 
ing the information needed to decide how the system 
should be perturbed to bring it to the target state. 
The theoretical tools required for obtaining this continu- 
ous estimate are those of continuous quantum measure- 
ment II H El El El El El IIS IM El 111 III III El. 

An estimate of the state of the quantum system (in our 
case the motional state of the atom), conditioned con- 
tinuously upon the results of measurement, is given by 
a stochastic master equation (SME). The concept of an 
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SME was first developed by Belavkin in [3 (as the quan- 
tum equivalent of the classical Kushner-Stratonovich 
equation), and a derivation of an SME in the language 
of quantum optics may be found in |56l | . In principle the 
observer may calculate this state estimate by integrat- 
ing the relevant SME using the measurement results, and 
from this determine the values of the inputs to the system 
to effect control. However, in practice the complexity of 
the SME will prevent the observer from calculating the 
state estimate in real time for systems where the dynam- 
ics are sufficiently fast, such as the atomic dynamics we 
consider here. As a result, a simplified state-estimation 
procedure is essential. 

The theory of feedback control for classicall y ob served 
linear systems is well developed [stI Issl. 59, 60]. This 
theory considers the optimality and robustness of feed- 
back algorithms for given resources, and in some special 
cases may be applied directly to the control of quantum 
systems H, . Discussions of the application of classical 
feedback control techniques to quantum systems may be 
found in Refs. Hljlllll. However, for nonlinear control 



be written in Ito form as [67 



problems [isl 161116^. such as that of the atomic motion 
we consider, no general results regarding the optimality 
of control algorithms exist as of yet. Here we will be 
concerned with presenting simple control algorithms for 
cooling the atoms, and demonstrating that along with 
an approximate estimation algorithm, a substantial in- 
tervention into the dynamics of the atom using real-time 
feedback control should be achievable experimentally in 
the near future. In doing so, we build upon previous re- 
sults by providing a detailed analysis of the cooling 
performance in this system, including its robustness to 
parameter variations and real-world limitations on com- 
putational speed and detection efficiency. Furthermore, 
we are able to anticipate several difficulties in an exper- 
imental implementation and suggest solutions to over- 
come them. 

Owing to the computationally expensive nature of sim- 
ulating an atom in an optical cavity, including the atomic 
motion in one dimension, supercomputers are invaluable 
for the task of evaluating the performance of the control 
algorithm over many realizations. This is the first time 
that full quantum simulations of such a CQED system, 
including the quantized, one-dimensional motion of the 
atom, have been performed. We are thus able to confirm 
the results of various approximations which have been 
used previously to provide a simplified analysis of this 
system [H H H El . 



II. DESCRIPTION OF THE SYSTEM 

The system we consider here is that of a two-level atom 
in a driven, single-mode optical cavity, where the out- 
put from the cavity is monitored using homodyne de- 
tection [sslls^ . Due to the continuous observation, the 
evolution of the system is described by a stochastic mas- 
ter equation (SME) for the density operator, which can 
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(1) 



dp = -i[{H/h), p]dt + -idt I iV(u)X>[CTe-'"^/'']p du 

J-Tik 

+KV[a]pdt + y/lfRn[a]p dW, 

where the Hamiltonian is 

^ = -^+5Cos(fca;)(cr^ae^'^* + H.c.)-£;(a + a^), (2) 
and we have defined the superoperators I?[c] and Ti.[c] by 



I?[c]p := cpc' — ^{c^cp + pc'c) 
T-i[c\p cp -I- pc^ — {c + c^)p 



(3) 



for an arbitrary operator c. In addition, the measured 
homodyne signal, given by 



dr{t) = riK{a + a))dt + ^Jr/RdW, 



(4) 



continuously provides information about the optical 
phase shift due to the atom-cavity system and hence, 
as we will see below, information about the atomic posi- 
tion. In the above equations, x is an operator represent- 
ing the atomic position along the length of the cavity, p 
is the conjugate momentum operator, a is the annihila- 
tion operator for the cavity mode, and a is the lower- 
ing operator for the atomic internal states. The energy 
loss rate of the cavity (or equivalently, the fuU-width-at- 
half-maximum cavity transmission) is k, and E gives the 
strength of the drivi ng laser, a nd is related to the laser 
power P hy E = ^JnPjJfujj^, where is the angular 
frequency of the laser light. We take the laser frequency 
to be resonant with the cavity mode. The laser angular 
wave number is given by = Wl/c, the decay rate of 
the atomic excited state is denoted by 7, g is the CQED 
coupling constant between the atomic internal states and 
the cavity mode, and A is the detuning between the atom 
and the cavity mode, given by A = — (cjl — ^^a), where 
cja is the atomic transition frequency. We will assume 
a positive (red) detuning, corresponding to a trapping 
optical potential. Also, N{u) is the probability that the 
atomic momentum recoil along the x-direction is u due 
to a spontaneous emission event, and dW is the differ- 
ential of a stochastic Wiener process. Note that we have 
written the stochastic master equation in the interaction 
picture where the free oscillation of the light and the 
dipole rotation of the atom are implicit. 

The density operator p(t) obtained as a solution to 
the SME condenses our knowledge of the initial condi- 
tion and the measurement record, allowing a prediction 
of the statistics of any future measurement on the system. 
As we show numerically in Section FlV CI the solution to 
the equation is, generically, insensitive to the initial con- 
ditions except at very early times. We therefore refer to 
the solution of the SME as the observer's best estimate 
or state of knowledge of the quantum system. 
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In order to perform numerical calculations, it is sen- 
sible to choose units for position, momentum and time 
that are relevant for the corresponding scales of the sys- 
tem. When the detuning of the atom from the cavity is 
large compared to the cavity decay rate and the CQED 
coupling constant, both the internal atomic states and 
the cavity mode may be adiabatically eliminated, and in 
this case the atom sees an effective potential, given by 



2 2 

Vcg{x) = —h^-^cos^{kx), 



(5) 



where a := 2E/k. Since this regime is the one in which 
we will be primarily interested, convenient units for x 
and p are given respectively by the width in position 
and momentum of the ground-state wave function in the 
harmonic approximation to one of the cosine wells. A 
convenient unit for time is set by the frequency of oscil- 
lation in the same harmonic potential. Thus, given the 
three scales of the problem. 



Wx 



P 



= agky^2h/\mA\ 



(6) 



we will use the scaled variables X := x/wx and P := 
p/wp, and we will measure time in units of 2tt/luiio- The 
SME for the density matrix may now be written as 



dp = -i[H, p]dt + jdt j_j^N {u)V{ae~'^^)p du 
+kV[a]p dt + ^n[a] dW, 



(7) 



where 

H = TTP^+g cos{kX){a^ae'^' + H.c) - E{a + a^) (8) 
and 



dr{t) = 7](a + a'')dt + ^ dW, 



(9) 



Thus, the atom moves in a sinusoidal potential, with pe- 
riod Tr/k. The average height of the potential wells is es- 
sentially proportional to (a^'a). It is therefore instructive 
to consider the cavity- mode dynamics, which determine 
the magnitude of the effective potential. Under the ef- 
fective Hamiltonian, the Heisenberg equations of motion 
for the relevant cavity mode operators are 



dta — iE 



— a + i— cos'^{kX)a + Vkain{t) 
2 A 



dt{a' a) ~ —ka) a — lEia — a)) 



(11) 



where ain(i) is the input quantum noise operator com- 
ing from the electromagnetic field outside the cavity, and 
satisfies [ain{t) , a\^{t + r)] = 5{t). We can now identify 
a regime in which k is large compared to g^/A. In this 
case the cavity mode damps quickly to a steady state, 
and we can approximate the solution by 



l-i2 



kA 



cos^{kX) 



la — a 



Ak 



cos^ikX), 



(12) 



(a 



2i 



There will be fluctuations about these values due to the 
noise with a spectrum with a width of the order of k. 
When K (or k) is sufficiently large, then (a^a) « a^, and 
the height of the potential will therefore be determined 
by the strength of the driving. As a consequence, one can 
manipulate the potential that the atom sees by changing 
the strength E of the driving. The measured homodyne 
signal is given in the adiabatic approximation by 



dfit) = -y/8^{cos'^ {kX))dt + y/TjdW, (13) 

where dr{t) := dr{t)/^/k; the position-information con- 
tent of this signal is more clear in this form. 



with dr{t) = dr{t)/^/k. The scaled rate constants 

K, -E, A,7, and g — \/ nA/{ak) are obtained from their 
absolute counterparts by dividing by oj-ao/i^Tr). The 
scaled wave number is naturally k = kwx- Note that in 
these scaled units we have scaled h out of the problem, 
so that [X, P] = i. 

While we will solve this equation numerically, it is nev- 
ertheless useful to examine the solutions obtained by adi- 
abatically eliminating both the atomic internal states and 
the cavity mode, which is possible when the detuning A 
is much larger than the coupling constant g so that the 
excited atomic state is nearly unpopulated during the 
evolution 66]. In this case the Hamiltonian part of the 
evolution is governed by the effective Hamiltonian 



~9' 



■■^cos\kX)a 



(10) 



III. COOLING THE ATOMIC MOTION USING 
FEEDBACK 

As mentioned in the introduction, in implementing 
feedback control, two considerations are paramount. The 
first is the ability to effectively track the evolution of the 
system in real time (state estimation) . The second is to 
provide an algorithm which determines how this informa- 
tion about the system should be used to alter the inputs 
to the system to effect control. We now consider these 
two questions in turn. 



A. State Estimation in Real Time 

The SME given in Eq. Q continuously provides the 
observer's best estimate of the state of the system from 
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the information provided by the measurement record r{t) 
[Eq. I0J]. Therefore, it is possible, in principle, to use the 
measurement record to integrate the full SME to obtain 
a continuous state estimate. However, the resources re- 
quired to do this in real time are prohibitive in practice. 
This is because the SME, being a partial stochastic dif- 
ferential equation, contains a large effective number of 
variables, and is almost impossible to integrate in real 
time, given that the time scale of the atomic motion is 
typically a few /iS. As a result, we must obtain an es- 
timation equation which contains only a small number 
of variables, but is nevertheless a good approximation to 
the full SME. 

To obtain a compact system of estimation equations, 
we first note that initially Gaussian states remain close 
to Gaussian for a fairly large number of oscillations in 
a sufficiently strong sinusoidal potential. We thus make 
a Gaussian approximation for the density matrix, which 
amounts to ignoring all cumulants higher than second 
order. In the present system, which effectively measures 
(cos^ (kX)), we expect this approximation to remain valid 
as long as the wave packet is much narrower than a single 
well {kAX <C 1), and the wave-packet momentum width 
is much larger than two photon recoils (AP ^ 2k). To 
further simplify the estimation equations, we consider 
the SME after adiabatic elimination of the cavity and 
internal atomic degrees of freedom, given by 



dp ^ -~i[Hcs, p] dt + 2TV[cos^ {kX)]pdt 
- y/2^n [cos^ (kX)]pdW, 



(14) 



which gives a three-parameter model {k, F, and 77). Here, 
we have defined the effective Hamiltonian 



where the maximum light shift 



A 



fc2 



(15) 



(16) 



is not an independent parameter in this static problem, 
but we will treat it as such when considering feedback 
via a time-dependent potential. We have also defined 
the effective measurement strength 



r := 



2a^ 

A2k 



(17) 



We further simplify the estimation model by making a 
Gaussian ansatz for the density operator. The resulting 
estimation equations contain merely 5 variables: the esti- 
mated mean position and momentum, {X)c and (P)o, the 
estimated variances, Vx and Vp, and the estimated sym- 
metrized covariance Co := {XP + PX)c/2 - (X)o(P)o- 
(Due to the information content of the measurement 
process, the two variances and the covariance are alge- 
braically independent, in contrast to a closed system.) 



X 



Using the equation of motion for the expectation value 
of an arbitrary operator A, 

d{A) = TT:[Adp] 

= -i{[A,H^s])dt + 2V{D[cos'^{kX)\A)dt (18) 
- ^/2^{n[cos^(kX)]A)dW, 

we can calculate the system of five Gaussian estimation 
equations, which are 

d(X)e = 27r{P)edt 

+y/»nrkV^ exp(-2fc2v;^) sm{2k{X)^)dWe 
d(P)c = -V;„axfcexp(-2pT/^)sin(2fc(X)e)dt 

+^/8^kCe exp{-2PV^) sin(2fc(X)e)dW; 
inCcdt 

-SrjrPV^^ exp(-4fc2T/^) sin^ {2k {X)c)dt 
+2y/&qrPV^^ exp{-2Pv^) cos{2k{X)c)dWe 
-4Kiaxfc^Cc exp(-2py^) cos(2fc(X)c)dt 
-FrP[l - exp{^8PV^) cos(4fc(X)e)]rft 
-S-nTk^C^ exp(-4fc2v;|) sin2(2fc(X)e)dt 

- 4(C2 + PV^) cxp(-2fc2y^)] 

X cos(2fc(X)o)dWc 
dCe = 2'KVlidt 

-2Vrm,^Pv^ exp(-2pF_^) cos{2k{X)^)dt 
-SrjVpV^Cc exp(-4fc2\/^) sin2(2fc(X)c)dt 
+2y/&fJrPV^C^ exp{-2Pv^) cos(2fc(X)e)dH/e- 

(19) 

Although we can eliminate T^nax from these estimator 
equations (since T^naxfc^ — tt), we keep the optical po- 
tential explicit here because both Knax and F will effec- 
tively become time-dependent quantities when we con- 
sider feedback control by modulating the amplitude of 
the optical potential. The reconstructed "Wiener incre- 
ment" dWc arises by making the Gaussian approximation 
to the photocurrent equation (|13|l : 



dV?. 



dWe 



df{t) 



exp(-2fcV_J) cos(2fc(X)e)]dt. 

(20) 

These estimation equations resemble stochastic differen- 
tial equations, but it is important to recognize that dWc 
does not have the same statistics as the usual Wiener in- 
crement dW . In particular, the reconstructed increment 
is the usual Wiener increment with an extra deterministic 
component. 



dW^ = dW + 



exp(-2fc2yj) cos(2fc(X)e) 
- (cos(2fcX)) 



dt, 
(21) 

which reflects the difference between the actual and es- 
timated quantum states. This is the route by which the 
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measurement information is incorporated into the esti- 
mator evolution. Thus, these equations cannot be treated 
by the usual high-order numerical methods for stochastic 
differential equations, which assume that the determin- 
istic and stochastic parts of the differential equation can 
be explicitly separated. 

Another important issue to note is that because the es- 
timator only has information about {cos^{kX)), it cannot 
tell us in which well the atom is located, and therefore, 
while feedback will allow us to cool an atom within a well, 
we cannot know within which well the atom is confined. 
Moreover, the estimator cannot tell us on what side of 
a well the atom is located. Note that this degeneracy 
is not a consequence of the Gaussian approximation, but 
also applies to an observer using the full SME for estima- 
tion. As a consequence, the estimator will show us the 
motion of the atom up to a phase factor of tt: the esti- 
mated state will either be in phase with the true motion, 
or completely out of phase, since these two motions are 
mirror images of one another. Thus, for the feedback al- 
gorithm to be effective for every run, it is important that 
it works for both cases without change. The algorithms 
that we describe in the next section have this property. 

It is worth noting that the insensitivity of the feed- 
back algorithm to the side of the well the atom is on 
has a potential advantage. In simulating the dynamics 
of the atom, we find that if the initial energy of the atom 
is large enough for it to move from one well to another, 
then as it does so the wave function tends to split across 
the interwell barrier, with part going back down into the 
original well, and part going down into an adjacent well. 
In this case the wave function is no longer Gaussian, but 
consists of two or more "humps." However, even though 
the two humps are on opposite sides of their respective 
wells, they have approximately the same downward mo- 
tion. Thus, even though the estimator is Gaussian, the 
fact that the humps are in different wells is immaterial to 
the measurement of (cos^(fcX)), and the estimator effec- 
tively tracks the evolution of an "equivalent" quantum 
state localized within a single well in the sense of gen- 
erating the same measurement record as the fragmented 
state. As a result, we can expect the Gaussian estimator 
to continue to work even though the atom is distributed 
across more than one well. 

In order to perform the continuous estimation 
(whether using the full SME or the above Gaussian ap- 
proximate estimator) one chooses an initial state for the 
estimator that reflects the initial ignorance of the true 
quantum state and integrates the estimator equations us- 
ing the measurement record as this record is obtained. In 
our case, when the atom is dropped into the cavity we 
have very little idea of where it is. It is therefore sensible 
to choose as the initial state for the estimator a Gaussian 
broad enough to cover one side of a well, and centered 
somewhere on that side (it is not important where pre- 
cisely). We only need to include one side of the well in 
our initial uncertainty, since as mentioned above, both 
sides are equivalent from the point of view of the estima- 



tor. In practice, though, it is sufficient to pick a state 
with a much smaller uncertainty product, centered on 
one side of a potential well. Even if this estimator state 
does not substantially overlap the actual quantum state, 
as we will see below, it will converge to the "true" state 
in a time of order l/F as the measurement information 
is incorporated. 

B. An Effective Cooling Algorithm 

As discussed in Section ^ it is possible to raise and 
lower the height of the potential by changing the input 
laser power. Thinking of the atom as a classical particle, 
one would expect to be able to cool the atom using the 
following simple algorithm: When the atom is moving to- 
ward the center of the well, it is being accelerated by the 
potential, and so we can lower the laser power to reduce 
this acceleration. On the other hand, when the atom is 
moving away from the center of the well, then the poten- 
tial decelerates the atom, and in this case we can increase 
the laser power so as to increase the deceleration. Writ- 
ing the optical field amplitude in the absence of feedback 
as E, we will denote the value of the field during times 
when it is "switched high" as (1 -t- £i)E and the value 
when it is "switched low" as (1 — e2)E. Thus the poten- 
tial height will have the respective values (1 +ei)'^V and 
(1 — £2)^^, where V is the unmodulated potential depth. 
In this way, we repeatedly switch the potential between 
a high value and a low value, which both slows the atom 
and moves it in towards center of the well, reducing the 
total energy and hence cooling the atom. This algorithm 
also has the nice property that all the information it re- 
quires is whether the atom is climbing or descending the 
well. Therefore, the question of whether the estimated 
state is in phase or completely out of phase with the 
true motion is immaterial, since in both cases the rele- 
vant information is known. However, it turns out that 
this algorithm is not very effective at cooling the atomic 
motion, because while it does cool the centroid, it simul- 
taneously feeds energy into the motion by squeezing the 
atomic wave function. A detailed discussion of this effect 
is given in Appendix A. The result is that this simple 
centroid-only cooling algorithm is not adequate to cool 
the atom to the ground state. To do so we must develop 
a more sophisticated algorithm which can reduce the mo- 
tional energy associated with the variances of the wave 
function. 

To proceed, we consider the change in the motional 
energy due to time dependence of the optical-potential 
amplitude: 

dt{E,s)fb = -(atKnax)(cos2(fcX)). (22) 

It is clear from this expression that the "bang-bang" feed- 
back strategy from the simple cooling algorithm, where 
the potential amplitude is switched cyclically and sud- 
denly between two values, maximizes the energy extrac- 
tion rate if we consider only algorithms with a cyclic 
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modulation of the potential within that range. In fact, 
we should switch the potential low when the value of 

— (cos^(fcX)) is maximized (i.e., the modulus is mini- 
mized), thus maximizing the energy extracted from the 
atom. Then we should switch the potential high when 

— (cos^(fcX)) is minimized, thus minimizing the energy 
transferred to the atom when completing the modulation 
cycle. 

When we must rely on the Gaussian estimator to de- 
termine when to switch the potential, we switch the po- 
tential low or high when the estimated quantity 

2/est := - exp(-2pF^) cos(2fc(X)e) (23) 

is maximized or minimized, respectively. Doing so in- 
volves an additional technical complication, however, 
since at any given instant, we must predict whether or 
not j/ost is an extremal value. To do this we dynamically 
fit a quadratic curve to a history of these values (typically 
the last several hundred), and then trigger the feedback 
transitions on the slope of the fitted curve. This pro- 
cedure is a convenient method for predicting the times 
when j/cst is maximized or minimized, and helps to reject 
the residual noise in the time evolution of this quantity. 
More formally, given a set of estimates yest,n correspond- 
ing to times t„, we implement the curve fit of the function 
ao + aix + a2X^ to the last q values of ?/ost,j- We will de- 
fer the discussion of the computational efficiency of this 
algorithm until Section [C2l 

It may seem strange that this algorithm makes use of 
past state information, when all the necessary informa- 
tion should in principle be contained in the current esti- 
mate of the state. Indeed, it is clear from the quantum 
Bellman equation Hit [T3 | that the optimal control al- 
gorithm is given by only looking forward in time given the 
current state. In practice, such an algorithm is difficult 
to realize, and the present algorithm is likely suboptimal 
but robust and effective. 

We also note that since we are triggering in this al- 
gorithm on the estimated value of — (cos^(A;X)), it seems 
that it could be much simpler to simply trigger on the de- 
tector photocurrent dr{t), which from Eq. H13|) we know 
is a direct measure of this quantity plus noise. In fact, 
we can use the same curve-fitting procedure on dr{t) 
to reject the noise on the signal. However, as we will 
see, this procedure does not work nearly as well as with 
the Gaussian estimator, because the signal df{t) is noise- 
dominated, and the Gaussian estimator acts as a nearly 
optimal filter for the useful information from this signal, 
in the same sense as a Kalman filter ,5^ ,6^ . 

C. Cooling Limits 

At this point, we can work out a simple theory of cool- 
ing for the improved cooling algorithm developed in the 
previous section. When we switch the optical potential 
low, the change in energy from Eq. (|22|l is given by 



where we have taken ei = £2 = £, so that the potential 
amplitude is switched from (1 -f e)^Vinax to (1 — e)^Vmax- 
When we switch high, the expression for the energy 
change is the same except for an overall minus sign. 

Now we consider these expectation values for a particle 
very close to the ground state, and we will first consider 
the case where all the excess motional energy is associ- 
ated with the motion of the wave-packet centroid. For 
the ground state itself, an equal part Eo/2 of the ground- 
state energy Eq is associated with each of the kinetic and 
potential parts of the effective Hamiltonian. If we define 
AE := (E'eff) — Eq as the small, additional motional en- 
ergy above the ground state, then as the displaced wave 
packet (coherent state) evolves in the nearly harmonic 
potential, the value T4iax([l — cos^(fcX)]) of the poten- 
tial part of the atomic energy oscillates between £^0/2 
and Eo/2 + AE. Then the best cooling that can be 
achieved for this state corresponds to two up and two 
down transitions per oscillation period (which is unity 
in our scaled units). Inserting these extremal values of 
— Vmax{cos^{kX)) into Eq. H24I) and the corresponding 
expression for upward transitions gives the cooling rate 
of 

A{E,s)tb = -8s{{E,s) - Eo) (25) 

per unit time. 

In practice, the cooling algorithm switches more often 
than in this estimate due to noise on the signal. However, 
in such cases we might still expect that Eq. H25|l is a 
reasonable estimate for the cooling rate. Even though the 
potential is switched more often, the switches are caused 
by noise and thus do not occur at optimal moments in 
time. Thus the amount of cooling per switching cycle is 
reduced, and we observe that the net cooling effect per 
unit time is near the optimal rate of Eq. H25|l , presumably 
because the cooling algorithm is not too sensitive to the 
noise. If we then assume that heating due to spontaneous 
emission is negligible compared to measurement heating, 
then the steady-state motional energy is obtained when 
the cooling rate in Eq. H25(l balances the heating rate 
which we derive in Appendix IbI fEo. IjBTp l. This gives 
the steady-state energy 

(i?eff)ss = (26) 

where /3 := rfc''/2£. 

The other extreme case we will consider is where the 
motional energy in excess of the ground-state energy is 
associated purely with squeezing of the wave packet and 
not with the centroid. In this case, Vinax(cos^(fcX)) os- 
cillates between the values {{Ecs) ± 1/ (-Ecff)^ — Eq)/2, 
giving the cooling rate of 

A{E,ff)n> - ~8e^E^^-E^ (27) 



A(i;eff)fb = -4£K„ax(cOs2(fcX)), (24) 



per vmit time. Then the condition that this cooling rate 
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balances the heating rate in Eq. (|B7|) impHes 



Eo 



(28) 



for the steady-state energy. But in the regime where 
these theories are vahd (i.e., (£'cff)ss — Eq <^ Eq or 13 <^ 
1), the centroid result of Eq. (|26() is always larger than 
the squeezing result of Eq. H28I) . Thus, we expect that the 
centroid limit [Eq. (|26|l ] is the more appropriate cooling 
limit. 

Eq. (|26|l predicts that the atomic motion will cool es- 
sentially to the ground state so long as /3 < 1/2. Thus, as 
various parameters are changed, there are "border" val- 
ues that divide the parameter space according to whether 
or not the atomic motion is cooled to the ground state. 
For example, the atom will cool to nearly the ground 
state so long as the "bang amplitude" e is larger than 
the border value 



£b = rfc^. 



(29) 



Otherwise, the steady-state energy should be substan- 
tially higher, and the atomic motion may not even cool 
at all. However, this simple theory is only valid at small 
energies, so such behavior should not necessarily be pre- 
dicted well by this theory. Additionally, this theory does 
not account for how well the Gaussian estimator tracks 
the true atomic state, and thus we would expect the ac- 
tual steady-state temperatures to be higher than pre- 
dicted by Eq. if^ . 

One final modification to this cooling theory is nec- 
essary, since as we discuss below, parity considerations 
show that half of the atoms can cool at best to the first 
excited band of the optical potential. To take this effect 
into account, we simply modify Eq. (|26|l by replacing Eq 
by the average of the ground and first excited band ener- 
gies, {Eq -\- Ei)/2. But in the harmonic approximation, 
Eq = TT and Ei = Stt, so the cooling limit becomes 



(Ecs) 



2-K 



(30) 



This modification is simple because the centroid argu- 
ment above only assumed that the atomic state is near 
the band with the lowest achievable energy. The odd- 
parity atoms have the first excited state as their lowest 
achievable state, and the above argument applies just as 
well to states near this effective ground state. Thus, the 
ensemble average amounts to simply averaging over the 
two possible steady-state energies. 



IV. SIMULATIONS: ADIABATIC 
APPROXIMATION 

A. Stochastic Schrodinger Equation 

Simulating the simplified dynamics in the adiabatic ap- 
proximation according to the SME in Eq. (|14|) is much 



less numerically intensive than simulating the full SME 
given in Eq. Q. Because of this, we will perform an ex- 
tensive analysis of the behavior of the feedback control 
in this section using simulations in the adiabatic approxi- 
mation. We will then verify that these simulations indeed 
provide a good description of the dynamics by perform- 
ing a limited number of simulations of the full SME in 
the following section. 

To perform our simulations, we first note that we may 
view any SME as being generated by averaging a stochas- 
tic Schrodinger equation (SSE), where the average is 
taken over all the signals that the observer fails to mea- 
sure 13 |63- The unobserved signals in our case are 
the spontaneous emission from the atom and the part of 
the cavity output that is not measured due to photode- 
tector inefficiency (77 < 1). We therefore replace nature 
with a formal omniscient observer, unknown to the ob- 
server who will effect control. The omniscient observer 
measures everything that the control observer does not, 
and also has access to the control observer's measure- 
ment record. From the omniscient observer's point of 
view, the evolution of the system is described by an SSE. 
We can simulate the system by integrating this SSE, and 
even though it does not provide us with the control ob- 
server's state of knowledge, which would require averag- 
ing over all possible measurement results of the omni- 
scient observer, it does generate instances of the control 
observer's appropriate measurement record. The control 
observer can then use this measurement record to obtain 
an estimate of the state of the system by integrating the 
Gaussian estimator derived above 0. Since the "cor- 
rect" state estimate, given by integrating the full SME, 
is not required, we can obtain a full simulation of the 
feedback control process by integrating the SSE, which 
requires far fewer resources than integrating the SME. 

The normalized Schrodinger equation corresponding to 
the master equation (|14|) that we integrate numerically 
is 



r 

-iH\ilj) dt - - cos^{2kX)\i}) dt 
-^{coB{2kX)) cos(2fcX)|V') dt 
{cos{2kX dt 



{cos{2kX)) ~ cos(2A:X) 



(31) 



<dW, 



which generates trajectories with the same measure as 
the SME. Since we are integrating the SSE in lieu of the 
SME, there are a few subtleties related to the measure- 
ment noise to be accounted for. In particular, the Wiener 
increment dW in the SSE 1)3 l|l is not equivalent to the 
Wiener increment in expression H13|) for the measurement 
record if 77 < 1, because the measured noise does not fully 
represent the "true" noise. Rather, the quantity y/rjdW 
in the photocurrent expression should be replaced by a 
measured noise dWrj, given in terms of the full noise of 
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the SSE by 



(32) 



where dVFaux is an auxihary Wiener increment, so that 



dW^ = rjdt, and a fraction 



dM^unm - (1 - 77) rfW^ - Vvi^ - V) dW^- 



(33) 



of the true noise is not measured by the control observer 
(i.e., dWrj+dWunm = dW). The photocurrent in Eq. ifT^ 
then becomes 

df{t) = -y^8^{cos^{kX))dt + dWr„ (34) 

and the estimated Wiener increment is given in terms of 
the measured noise as 



dW,, 



dW„ 



2f]T exp(-2fcVx)cos(2fc(X)o 
- (cos(2fcX)) 



dt, 
(35) 

which still reflects the difference between the actual and 
estimated quantum states in the same way as before. 



B. Details of the Simulations 

For the simulations, we choose a "canonical" set of 
experimentally realistic parameters that will put us in 
the regime where adiabatic elimination of the cavity and 
internal atomic states should be valid. Using the Ce- 
sium D2 line as the atomic transition, we have m — 
2.21 X 10-25 kg and Wq/Stt = 351.7 THz; we assume that 
the cavity subtends a small solid angle and thus that 
the spontaneous emission rate is given approximately by 
the free-space value of 7/27r = 5.2 MHz. For the cav- 
ity parameters, we choose values similar to those used in 
recent CQED experiments , yielding an energy 

decay rate of k/2tt — 40 MHz, an atom-field coupling 
constant of g/2TT = 120 MHz, and a mean intracavity 
photon number of a = 1. However, we will consider a 
much larger detuning of A/27r = 4 GHz to the red of the 
atomic resonance than was used in the experiments, in 
order to work in a dispersive regime where the adiabatic 
approximation should work well. The scaled parameters 
corresponding to these physical values are F = 23.6 and 
k = 0.155. The scaled depth of the optical potential is 
Knax = 7r/fc2 = 131, which corresponds to an atomic 
speed of 14.7 cm s"^ and is of sufficient depth that the 
lowest two band energies Eq — 3.12 and Ei — 9.33 are 
close to the corresponding values in the harmonic ap- 
proximation of TT and 37r, respectively (with 27 trapped 
bands in the optical potential and a ground-band width 
of 1.7 X 10-^2^ go that tunneling transitions between wells 
are highly suppressed). Also, for the canonical parameter 
set we choose a feedback amplitude e = 0.1 (easily ac- 
complished experimentally) and the idealized detection 
efficiency of 77 = 1. The goal of our analysis will be 



to understand the dynamics of the system under these 
ideal conditions and then to understand how variations 
in these parameters influence the dynamics. 

The optical potential in the simulations spanned 24 
wells, corresponding to the cavity lengths in the Cal- 
tech experiments [T^. To crudely model the sticking 
behavior of an atom when colliding with a mirror, we 
implemented absorbing boundary conditions, where the 
absorption ramped smoothly on in the last half of each of 
the boundary potential wells. We explicitly conditioned 
the simulations on the fact that they were not lost, mim- 
icking the experimental situation where a run with a lost 
atom would simply be repeated. However, given the fact 
that the cavity is relatively long compared to the opti- 
cal potential period, these details do not substantially 
influence the results of the calculation. 

To understand the typical dynamics of this problem, 
we consider ensemble averages. We consider the experi- 
mental situation where atoms dropped one at a time have 
substantial motional energy after loading into the opti- 
cal potential, but only atoms that are well trapped are 
candidates for further cooling. Thus, we assume that all 
atoms in the simulation ensemble have a centroid energy 
of 84.2, which corresponds to a stationary atom at a dis- 
tance of 6 away from the bottom of a well (where the 
well period is n/k ~ 20.3), and we distribute the initial 
locations of the atoms uniformly between the bottom of 
the well and this maximum distance. The atoms in the 
simulation begin in a coherent state with T4; = = 1/2 
in the center well of the cavity. We do not expect that the 
steady-state results are sensitive to the details of these 
initial conditions. Except where noted otherwise we plot 
averages over 128 trajectories; in plots where we inspect 
the variation of ensemble averages as a function of a pa- 
rameter, each ensemble average is computed with respect 
to an independent set of random- number seeds. 

We obtained numerical solutions to the SSE (|31|l using 
the order 1.5 (strong), implicit, stochastic Runge-Kutta 
(SRK) algorithm in Ref. JTlj as part of an operator- 
splitting method. To evolve the system over a large step 
of At = 0.0005, we first used the SRK evolver to evolve 
the wave function according to the SSE without the ki- 
netic 7rP2 term in substeps of size At/ 8 over an interval 
At/2. Next, we applied the operator exp{—iTTP^At) af- 
ter a Fourier transform to momentum space. Finally, af- 
ter transforming back, we again evolved according to the 
SRK algorithm by another interval At/2 (in four steps) 
to complete the full At step. This operator splitting pre- 
serves the order 1.5 (temporal) convergence of the SRK 
approximation. The spatial grid contained 2048 points. 
The space and time discretizations were found to give ad- 
equate convergence when computed in 32-bit precision for 
most of the trajectories, although a minority (typically 
only a few percent) displayed some sensitivity to these 
numerical parameters due to their proximity to unsta- 
ble hyperbolic points in phase space. However, for the 
cooling simulations, the steady-state energies were well 
converged and not affected by this sensitivity, since the 
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FIG. 1: (Color online) Example of locking behavior of the 
Gaussian estimator for centroids and variances. The heavy 
lines (a) are the expectation values for the "true" wave packet, 
while the light lines (b) are the corresponding Gaussian esti- 
mator quantities. The improved feedback cooling algorithm 
was switched on at t = 2. The bottom graph shows the corre- 
sponding photodetector signal increments Af, plotted for the 
time increments of Ate — 0.0005 used by the estimator. 



atoms only saw the harmonic portions of the potential 
wells. 

The Gaussian estimator begins in each case with the 
impure-state initial conditions (^)o = 6, {P)c — 0, 
Vx — Vp = l/\/2, and Co = 0. As we argued above 
in Section llll Al this choice of the initial estimate is arbi- 
trary but sufficient for our purposes, and we will see in the 
next section that the coohng performance is insensitive 
to any "reasonable" initial estimate. The estimator is up- 
dated with a time step of At^ = At using the simple Eu- 
ler method [tJ (both to keep the calculations simple and 
because, as we discussed in Section llll Al dWc is not the 



usual Wiener increment and thus violates assumptions 
used in constructing higher-order integrators). Some ex- 
tra measures are needed to guard against instabilities in 
the estimator evolution, which are due to two possible 
reasons. First, the low-order numerical integration of 
the stochastic estimator equations with a relatively large 
time step tends to be unstable, especially when the occa- 
sional large fluctuation occurs. Second, while Gaussian 
approximations to Hamiltonian evolution (or even evo- 
lution under an unconditioned Lindblad master equa- 
tion) can be thought of as resulting from a variational 
principle, and hence inherently stable, the Gaussian ap- 
proximation to the conditioned SME cannot, and thus 
does not inherit these strong stability properties. We ex- 
plicitly detect the four unphysical conditions: negative 
position variance (VJ < 0), negative momentum vari- 
ance (Vp < 0), phase-space area below the Heisenberg 
limit (Ae := [V^V^ - C^]!/^ < 1/2), and large position 
variance (V^ ^ 1, indicating the atom is not localized 
within the cavity). If any of these conditions are de- 
tected, the variances are reset according to the initial 
values Vx = Vp = Ce = 0, while the means {X)c 

and {P)e are not modified; the measurement information 
will quickly restore the correct estimator values. For a 
typical cooling simulation with the canonical parameters, 
the estimator must be reset in about 95% of the trajecto- 
ries, with trajectories on average requiring 3 resets due to 
the area condition (the other conditions producing neg- 
ligible numbers of resets in comparison). Also, to help 
alleviate the number of resets caused by this condition, 
we instead implement the condition that resets only oc- 
cur when Aq < 1/4. This situation may seem curious, as 
the Gaussian area evolves according to 

-SriTk^V^A^ exp{-Ak^V^)sm^(2k{X)^)dt 
-8r]Tk*V^'^[l - Ak^V^ eM-'^k^Vx)] 

X cxp{~2Pv^)coa^{2k{X)^)dt 
-^/2qrPv^[l - A{PV^ + Ac)exp(-2py^)] 

X cos^{2k{X)^)dWc, 

(36) 

whence it follows after some examination that A^ > 1/2 
for all time if this is true initially. Thus, these resets 
are numerical artifacts of evolving the estimator using 
a finite time step and a low-order numerical method; in 
the simulations, reducing the time step Ata by a factor 
of 10 reduces the average number of resets due to the 
area condition to 1.1 (with 0.8 resets on average due to 
the Vx < condition). However, we continue to use the 
larger step size, which is more reasonable to implement 
in practice given technological constraints on the speed 
of calculations for the forseeable future. For somewhat 
smaller measurement efficiencies, the estimator requires 
fewer resets (e.g., about 12% of the trajectories for rj = 
0.8, and only 1.7 resets on average per trajectory), since 
the steady-state area is larger (corresponding to a state 
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FIG. 2: (Color online) Evolution of the mean energy {Hcs} 
(relative to the minimum potential energy), comparing the 
performance of different cooling algorithms. Heavy solid line 
(a): improved cooling algorithm. Solid line (b): centroid-only 
cooling algorithm. Dashed line (c): no cooling algorithm. 
Each line is an average of 128 trajectories; parameters are for 
the canonical set. 

with lower purity, as Tr[p^] = l/{2Ae) [t^, where pc is 
the estimated density operator). 

Because the estimator locks on quickly with the canon- 
ical parameters, the feedback cooling algorithm is turned 
on at time t = 2. For the improved cooling algorithm, 
we find that fitting to the last 300 observations (spaced 
apart in time by Ate) is optimal. 

C. Estimation and Cooling Dynamics 

Now we will evaluate the performance of the two main 
components of the quantum feedback cooling system, the 
estimator and the cooling algorithm. Henceforth, we will 
only study the performance of the improved cooling al- 
gorithm, except where noted otherwise. Fig. ^ shows a 
typical trajectory undergoing cooling and compares some 
of the "true" moments computed from the atomic wave 
function (i.e., those computed with respect to the result 
of integrating the SME of the omniscient observer) with 
the moments from the estimator. The centroids and vari- 
ances of the estimator lock on by the time that the cooling 
algorithm is switched on (t = 2), and the estimator pro- 
vides a reasonably faithful representation of the atomic 
dynamics. It is especially important that the variances 
track the atomic wave packet accurately in order to make 
the improved cooling work effectively. Fig. ^ also shows 
the measurement record increments 

pt + AtG 

Af(i) J df{t) (37) 

corresponding to the time increments A^q — 0.0005 used 
in the Gaussian-estimator integration. The measurement 
record appears to be noise-dominated, but nevertheless 
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FIG. 3: (Color online) Evolution of the mean energy {H^g}, 
comparing the performance of the improved cooling algorithm 
with different initial (Gaussian) state estimates. Heavy solid 
line (a): estimate matches the canonical parameter set. Solid 
line (b): initial estimate is much more impure and substan- 
tially overlaps {x,p) = (0,0). Dashed line (c): same as solid 
line, but centered on {x,p) = (0,0). 



contains information about the quantum state in the 
time-dependent offset level of the noise. The action of 
the estimator is such that it "demultiplexes" the infor- 
mation about the motions of the centroids and variances, 
even though they are encoded in the same frequency 
range in the photocurrent signal. One can get a sense 
for how this works from the form of Eqs. (|19l) . For exam- 
ple, the measurement (dWc) terms in the d{X)c and dV^ 
equations are weighted by sin(2fc(X)e) and cos{2k{X)e), 
respectively. Thus, the incoming measurement informa- 
tion influences {X)^ or V^, depending on the state of the 
estimator: if the estimator is centered in a well, varia- 
tions in {cos^ kX) (i.e., the photocurrent) are presumed 
to be caused by variations in Vx, while such variations 
are attributed to {X) if the estimator is centered on the 
side of a well. The estimator can also gain information 
from the absolute value of the measurement record — as 
we discussed in Section IIII CI energy in only the vari- 
ance degree of freedom can produce transient values of 
14iax([l — cos^(fcX)]) below Eq/2, whereas energy in only 
the centroid cannot. This subtle unraveling of the mea- 
surement information into the estimated state is auto- 
matic in the measurement formulation we have used here. 

Fig.[21compares the ensemble-averaged evolution of the 
energy (-Eeff) in the case where the atoms are not cooled 
to the cases where the atoms are cooled by the centroid- 
only and improved cooling algorithms. When the atoms 
are not actively cooled, the energy simply increases lin- 
early in time. The initial heating rate of dt{Ecff) = 1.78 
predicted by Eq. ljB6p for measurement-backaction is sub- 
stantially smaller than the initial simulated heating rate 
of 2.1(1). We can understand this discrepancy, since 
the initial conditions in the simulation are uniformly dis- 
tributed in phase space along a constant-energy surface, 
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FIG. 4: (Color online) Evolution of the mean energy (_ffcff)? 
comparing the performance of the improved cooling algorithm 
(heavy solid line, a) with that of cooling based directly on the 
homodyne signal (solid line, b). Heavy solid line: improved 
cooling algorithm. Each line is an average of 128 trajectories; 
parameters are for the canonical set. 



weighting more heavily the gradients of the potential 
where heating is maximized; Eq. ljB6p underestimates the 
heating rate by assuming a uniform distribution in con- 
figuration space. 

The simple, centroid-based cooling algorithm cools the 
atoms rapidly for short times, but around t — 10, the 
cooling stops and the atoms begin to heat back up, ap- 
parently without bound. Clearly, we expect some heating 
from the squeezing effects discussed above. Furthermore, 
we expect that when the centroid cools and the variance 
becomes much larger than the mean-square amplitude 
of centroid motion, the centroid evolution will be dom- 
inated by measurement (projection) noise, and thus the 
feedback signal will also be determined mostly by random 
noise from the measurement. In this state, the cooling 
algorithm is ineffective and leads to heating. This pro- 
cess runs away, since there is no mechanism to cool the 
energy associated with the variances (squeezing). 

On the other hand, the improved cooling algorithm 
also rapidly cools the atoms and the atomic energies set- 
tle to a well-defined steady state, which already demon- 
strates the potential utility of this cooling algorithm for 
long-term storage of the atoms. For the moment, we will 
defer the question of how closely the atoms are cooling 
to the lowest energy band of the optical potential until 
we take a closer look at the atomic dynamics in the next 
section. 

The cooling behavior that we see is also insensitive to 
the exact choice of the initial (Gaussian) state estimate, 
so long as we have made a reasonable choice. To illus- 
trate this statement, Fig.|21shows the ensemble-averaged 
cooling dynamics for two initial state estimates in addi- 
tion to the canonical estimate that we use for the rest of 
the simulations. In the first modified initial estimate, we 
have selected (X)e = 0.1, (P)o = 0,V^ = 1, = 1, and 



Ce = for a much larger uncertainty product (lower pu- 
rity Tr[/9^]) and a substantial overlap with the other side 
of the potential well. We see that although the short- 
time cooling performance is slightly worse, the long-time 
cooling performance is the same as for the original initial 
estimator. This implies that the tracking behavior is ini- 
tially worse for this initial estimator, but at sufficiently 
long times the estimator "forgets" its initial state and 
does a better job of tracking. The other initial estima- 
tor is the same as the previous case, but with {X)c = 0, 
so that the estimator is centered on the origin in phase 
space. Examination of the estimator equations H19|) re- 
veals that in this case the centroids {X)c and {P)c will 
remain zero for all time, effectively "freezing out" the 
estimator's centroid degree of freedom. The cooling be- 
havior here is much worse than in the other two cases, 
demonstrating both that the estimator's variance degree 
of freedom is insufficient to mimic the dynamics of the 
true state and that it is important to exercise some cau- 
tion when selecting an initial state estimate. 

Finally, we return to the question raised in Sec- 
tion nil Bl regarding the necessity of using the Gaussian 
estimator at all, in light of the fact that the measure- 
ment record itself provides a direct^ albeit very noisy, 
measurement of the quantity (cos^ kX) used to deter- 
mine the control switching times in the improved cooling 
algorithm. Fig.^jaddress this question by comparing the 
energy evolution under the improved cooling algorithm 
using the Gaussian estimator with the same algorithm 
but based directly on the measurement record. In both 
cases, the quadratic curve was fitted on-the-fly to the last 
300 points in the measurement record to reduce noise ef- 
fects. We found this number of points to be optimal for 
the canonical parameter set, and this optimal time inter- 
val over which the data are fit is tied directly to the time 
scales of the atomic motion. Although the direct-cooling 
case cools rapidly and settles to a lower temperature, it is 
clear that the cooling performance is greatly improved by 
the Gaussian estimator. The estimator acts as a nearly 
optimal noise filter, efficiently extracting the relevant in- 
formation from the noisy measurement record. 

D. Parity Dynamics 

We now consider the evolution of the parity of the con- 
ditioned atomic state due to the measurement. For this 
discussion, we will focus on the expectation values of the 
usual parity operator V (where Vip{x) = ip(—x)). For 
numerical purposes, though, we will instead use the re- 
duced parity operator, V' := TZ^VTZ, where the reduction 
operator TZ is defined by 

R,W:= f;*(.-f),«[-i^), (38) 

J — — OG \ / L / 

which allows us to detect the parity of the atomic state 
with respect to any potential well in the cavity, rather 
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FIG. 5: (Color online) Histograms of the reduced parity {V') 
for an ensemble of 2048 trajectories undergoing cooling, show- 
ing how parities evolve towards the extreme states of pure 
parity at late times. 



than just the center well. We will simply use the notation 
'T" for both operators in the discussion, though, as the 
discussion holds for either operator. The initial condi- 
tions for the simulations, which correspond to Gaussian 
wave packets displaced in phase space from the origin in 
phase space, represent equal mixtures of odd and even 
parity {{V) = 0) to a good approximation. As time 
evolves, however, the parity of the conditioned atomic 
state evolves (in the adiabatic approximation) according 
to 



d{V) 



V8r/r {Vcos{2kX)) - {V){co8{2kX)) dW. 

(39) 

This effect is due only to the nonlinearity of the measure- 
ment term in the master equation H14|l : the parity under 
unconditioned evolution is invariant. 

The reduced parity evolution for a simulated ensem- 
ble of trajectories undergoing measurement and cooling 
is shown in Fig. |S1 It is clear that (V) evolves towards 
the extreme values ±1 of parity for the simulated tra- 
jectories. This "parity purification" behavior is not im- 
mediately clear from the evolution equation l|39|l . which 
states that {V) evolves according to a diffusion process. 
The key point is that this diffusion process is nonstation- 
ary, since the diffusion rate depends on the parity, and 
indeed vanishes for states with pure parity. To see the 
final steady state directly, we consider the projectors V+ 



and V- for even and odd parity, respectively, defined by 

l±V 



V± 



(40) 



Then the product of the expectation values of these pro- 
jectors evolves according to 

d[{V+){V-)] = -^r^T{{V+){V-)pfdt, (41) 

where the parity difference p is given for states with 
{V±) ^ by 



P 



{V+ cos2 kX) {V- cos^ kX) 



{V+) 



{V-) 



(42) 



The magnitude of p represents how well the measure- 
ment of cos^ kX can resolve the parity of the atomic 
state, and it is nonzero for generic states. Thus, the 
quantity {V+){V^) damps monotonically to zero, which 
is only possible if the parity of the state becomes either 
pure even or pure odd. Given the initial conditions in 
the simulations, and the fact that the diffusion process 
in Eq. H39|l is unbiased, we can expect that half of the tra- 
jectories will become even and half will become odd. As 
we mentioned above, this argument also applies to both 
the standard and reduced parity operators. However, the 
damping rate in Eq. (|41|l is greatly reduced for the stan- 
dard parity operator when considering states localized in 
potential wells not centered on X — 0. 

This purification effect obviously has consequences for 
cooling the atomic motion. Even if the cooling algorithm 
can, in principle, cool the atomic motion faster than it is 
heated by the measurement, the cooling algorithm cannot 
affect the parity of the atomic state, and thus odd-parity 
states will not be cooled to the ground state of the optical 
potential. At best, then, the cooling algorithm can cool 
the atom to the ground energy band with probability 1/2 
and the first excited energy band with probability 1/2. 
This cooling behavior is illustrated in Fig. El where the 
evolutions of the band populations are plotted for the two 
cases where cooling is based on the Gaussian estimator 
and directly on the true atomic wave function using the 
improved algorithm in both cases. We can clearly see 
that in steady state the lowest and first excited bands 
are about equally populated, with these bands account- 
ing for 94% of the population in the Gaussian-estimator 
case and 98% of the population in the perfect-knowledge 
case. Thus, this parity effect is the dominant limit to 
the atomic motional energy for the ensemble. However, 
this situation is essentially as good as cooling fully to 
the ground band, because the measurement can distin- 
guish between the two possible outcomes (i.e., by resolv- 
ing the two different possible steady-state potential en- 
ergies) . Experimentally, if the odd-parity outcome is de- 
tected, the atom can be transferred to the ground band 
by a coherent process (e.g., resonant modulation of the 
optical potential by an external field or two-photon, stim- 
ulated Raman transitions), or the state can simply be 
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FIG. 6: (Color online) Evolution of the band populations for 
the lowest two (0 and 1) and next two (2 and 3) energy bands 
of the optical lattice, averaged over 128 trajectories. Top: 
cooling based on the Gaussian estimator. Bottom: cooling 
based on perfect knowledge of the actual wave function. 



rejected, in hopes of obtaining the even-parity outcome 
on the next trial. 

Finally, we comment on one odd aspect of the 
intermediate-time histogram (t = 40) in Fig. [S] which 
shows a marked asymmetry, even though the median 
value of {V) is zero. Because we can write Eq. (|39|l in 
the form 



d{V) 



(43) 



all of the asymmetry is due to the p factor. Since this 
diffusion process mimics a quantum nondemolition mea- 
surement of the parity, we can interpret p as an effec- 
tive measurement strength that is larger for even-parity 
states than for odd states. This asymmetry is reason- 
able since the measurement and cooling processes pro- 
duce even-parity states that are more localized than the 
odd ones. Thus, the even-parity states, being tightly lo- 
calized at the field antinodes, produce more extreme val- 
ues of the measurement record and are hence more easily 
distinguished from mixed-parity states. Again, however, 
this effectively asymmetric parity-measurement process 
does not affect the long-time outcome that either parity 
will be selected with equal probability. 
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FIG. 7: (Color online) Final energies, measured over the inter- 
val from t = 90 to 100, as the switching amplitude e varies; 
other parameters match the canonical set. Circles: cooling 
based on the Gaussian estimator. Squares: cooling based on 
perfect knowledge of the actual wave function. Dashed line: 
simple cooling theory, Eq. 13011 . Inset: magnified view of the 
same data. Error bars reflect standard errors from averages 
over 128 trajectories. 
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FIG. 8: (Color online) Final energies, measured over the in- 
terval from f = 90 to 100, as the unmodulated optical po- 
tential depth Knax ~ Tv/k^ varies; other parameters match 
the canonical set. Circles: cooling based on the Gaussian 
estimator. Squares: cooling based on perfect knowledge of 
the actual wave function. Dashed line: simple cooling theory, 
Eq. 1301 . Error bars reflect standard errors from averages over 
128 trajectories. 



E. Tests of Cooling-Limit Theories 

Now we can examine the dependence of the steady- 
state cooling energies on the system parameters and test 
the validity of the cooling- limit theories in Section IlII CI 
Fig. [3 shows the dependence of the simulated steady- 
state temperatures on the control-switching amplitude e 
in both the Gaussian-estimator and perfect-knowledge 
cases. The simple cooling prediction of Eq. H3()(l is 
also shown here for comparison. The agreement be- 
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tween the theory and the perfect-knowledge simula- 
tions is much better than between the theory and the 
Gaussian-estimator simulations, which is sensible be- 
cause the simple theory does not attempt to account for 
the imperfections in how the estimator tracks the true 
state. While the (perfect-knowledge) steady-state ener- 
gies consistently lie slightly below the theoretical predic- 
tion, the agreement is overall quite good, and the theory 
correctly predicts a sharp transition from efficient to in- 
efficient cooling near the transition border — 0.014 
(corresponding to the point /3 ~ 1/2 where the predicted 
steady-state energy is twice the smallest predicted value) . 
Simply put, when the control transitions of the optical 
potential are sufficiently weak, the cooling algorithm does 
not cool the atoms at a rate sufficient to counteract the 
measurement heating (the system is no longer "control- 
lable"), and so the atoms do not equilibrate at a very 
low temperature. The cooling and transition behaviors 
are very similar in the Gaussian-estimator simulations, 
but these energies are overall slightly higher, again prob- 
ably due to the inability of the Gaussian estimator to 
perfectly track the atomic state, even with rj = 1. 

Similar behavior appears as the (unmodulated) optical 
potential depth VJnax = 7r/fc^ varies, as shown in Fig. |S1 
For sufficiently large potential depth, the temperature 
is essentially constant, but if the potential depth is too 
weak, we again lose controllability and the cooling algo- 
rithm fails to be effective. There is a much more sub- 
stantial difference between the Gaussian-estimator simu- 
lations and the simple cooling theory for small Vmax com- 
pared to Fig. [71 This discrepancy is most likely due to 
the reduced confinement of the wave packet, which pro- 
duces stronger anharmonic effects and thus a breakdown 
of the Gaussian approximation. The simulation data are 
still in good agreement with the simple theory. 

The dependence of the steady-state energy on the mea- 
surement strength F is shown in Fig. O In this case, 
the theory again matches the perfect-knowledge situa- 
tion quite well. The cooling is best for moderately small 
values of F, and the simple theory predicts a transi- 
tion to an uncooled state if F exceeds the border value 
Fb = e/k'^ = 173, since F controls the rate of heating 
in the system. However, due to the functional form of 
Eq. I|30|) , the "transition" here is much more gradual than 
the transition observed when varying e. The Gaussian- 
estimator simulations again result in values that are con- 
sistently slightly above the perfect-knowledge energies. 
The difference is very pronounced for small F, where 
we might expect that the estimator is not gaining suf- 
ficient information from the measurement and cannot 
track the true state well enough to produce effective cool- 
ing (i.e., the system is no longer "observable"). The 
perfect-knowledge energy data also show a slight upturn 
for small F that is not predicted by the simple theory. 
This departure indicates the direct importance of the lo- 
calization produced by the measurement in the cooling 
process: a localized phase-space distribution is responsi- 
ble for larger fluctuations in (cos^ kX) than a delocalized 
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FIG. 9: (Color online) Final energies, measured over the 
interval from t = 90 to 100, as the effective measurement 
strength F varies; other parameters match the canonical set. 
Circles: cooling based on the Gaussian estimator. Squares: 
cooling based on perfect knowledge of the actual wave func- 
tion. Dashed line: simple cooling theory, Eq. (1301 . Error bars 
reflect standard errors from averages over 128 trajectories. 
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FIG. 10: (Color online) Final energies as in Fig. |7| but with 
initial position and momentum centroids (for both the atomic 
wave packet and Gaussian estimator) at zero. The dashed line 
now refers to the variance-dominated, even-parity theory of 
Eq. (HHJ- 



distribution with the same energy, and from the discus- 
sion in Section IlII Bl the cooling rate is directly propor- 
tional to the magnitude of these fluctuations. Thus, a 
sufficiently large value of F must be chosen to maintain 
a localized atomic distribution and hence good cooling. 

Thus far, we have focused on testing the theory of 
Eq. (|30() . which is the most relevant theory for a physi- 
cal implementation of the cavity QED feedback system. 
However, we can also test the variance-only theory of 
Eq. (|28|) in the simulations, thereby giving indirect sup- 
port to the more physically relevant theory. We can do 
this by changing the initial condition of the atomic wave 
packets to be a coherent (Gaussian) state centered on 
one of the potential wells with zero centroid momentum. 
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From the form of the Gaussian estimator equations H28() 
we see that the wave-packet centroid will remain fixed 
at zero position and momentum, and thus all the en- 
ergy will be associated with the wave-packet variances 
(in the harmonic approximation). The results of sim- 
ulations of this type, where we also modify the initial 
condition of the Gaussian estimator to have the same 
centroid, are plotted in Fig.^| The theory of Eq. (|28|l is 
also shown; note that we do not need to replace the Eq 
in this expression with the average of Eq and since 
the initial conditions have exactly even parity, which is 
then time- invariant. The theory is again in good agree- 
ment with the perfect-knowledge data, both in the lower 
temperatures for large e and the sharper transition and 
lower transition border of £b = 0.0079 (corresponding to 
(3 ~ \/3/2, again when the predicted energy is twice the 
minimum predicted value). 

Finally, a note about the steady-state temperatures is 
in order for the cases where the Gaussian estimator is 
used for cooling. The ensemble-averaged, late-time en- 
ergies in Figs. I7I1UI for Gaussian-estimator cooling are 
slightly but consistently higher than the corresponding 
perfect-knowledge energies, suggesting that the cooling 
is somewhat less effective when the Gaussian estimator is 
used. However, these averages mask complicated dynam- 
ics that are illustrated in Fig. Trajectories that have 
even parity at late times settle into a quiet equilibrium 
with nearly all the population in the ground energy band. 
On the other hand, trajectories that settle to odd parity 
have a more complicated equilibrium, where the atom is 
mostly in the first excited band but the evolution is punc- 
tuated by "bursts" of heating. This is an indication that 
in steady state, the Gaussian estimator is, unsurprisingly, 
a much better approximation for a wave packet close to 
the ground state than a wave packet close to the first ex- 
cited state. As seen in Fig. ^2 these heating episodes cor- 
respond to population being transferred to higher energy 
states of odd parity, a process distinct from the parity 
drift which, after the initial cooling, transfers population 
between the ground and first excited states. Again, if the 
goal of the cooling process is to prepare an atom in the 
ground state, the rate of success is even better than sug- 
gested by the ensemble-averaged energies, since cooling 
in the even-parity cases is essentially as good as if one 
had perfect knowledge of the actual wave function. 



V. SIMULATIONS: FULL ATOM-CAVITY 
DYNAMICS 

To verify the validity of the results of this paper, it 
is important to check the adiabatic approximation that 
leads to the reduced master equation l(Tl|l . We thus car- 
ried out simulations of the full, coupled atom-cavity dy- 
namics described by Eq. (TJl. 

The simulation procedure is the same as before, ex- 
cept that we used the normalized Schrodinger equation 
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FIG. 11: (Color online) Evolution of the ground and first- 
excited band populations for three trajectories, illustrating 
typical behavior. Trajectories that purify to even parity (top) 
tend to settle to a quiet equilibrium in the ground state. Tra- 
jectories purifying to odd parity (middle) are characterized 
by dynamical equilibria, where periods in the first excited 
state are interrupted by episodes of higher energy. Trajecto- 
ries that take longer to purify in parity (bottom) cool rapidly 
to the lowest two states, after which population is transferred 
between these two states as the parity diffuses. Parameters 
here correspond to the canonical set. 

corresponding to the master equation {Tj), 

d\i}) = ~iH\ip) dt - '^a)a\il}) dt + '^{a + a))a\ii) dt 
(a + at)2|^) dt+^- {{a^a) - a) |^) dt 

+\/ka\ij) dW - ^(a + a^)|V') dW 

(44) 

We used the canonical parameter set described before. 
This set leads to scaled parameter values of 7 = 190 
and K = 1460 in addition to those used in the adiabatic- 
approximation simulations as described in Section llV Bl 
The detuning A/27r — 4 GHz from atomic resonance 
implies a wide separation of timescales in the prob- 
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FIG. 12: (Color online) Ensemble energy evolution comparing 
the effect of the spatial grid size (cavity mirror separation) on 
the atomic evolution. Heavy solid line (a): canonical parame- 
ters (adiabatic approximation) spanning 24 optical potential 
wells. Light solid line (b): 6 potential wells. Each curve rep- 
resents an average over 128 trajectories. 
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FIG. 13: (Color online) Ensemble energy evolution compar- 
ing evolution with and without invoking the adiabatic ap- 
proximation. Heavy solid line (a): adiabatic approximation. 
Light solid line (b): full atom-cavity dynamics at 4 GHz de- 
tuning from atomic resonance. Dashed line: full atom-cavity 
dynamics at 2 GHz detuning from atomic resonance. Each 
curve represents an average over 128 trajectories. 



lem. Thus we again evolved the system over a large 
step of At — 0.0005 but in much smaller substeps of 
Ai/4096, applying the kinetic evolution operator after 
half of the substeps. The spatial grid contained 512 
points, spanning only 6 wells to make the computation 
more tractable. The impact of the smaller grid on the 
cooling was to slightly accelerate the cooling process, as 
illustrated in Fig.^J This effect due to an effective evap- 
oration of hot atoms, since we treat the mirrors as ab- 
sorbing boundary condition as discussed above. We also 
used 7 cavity states and 2 internal atomic states. Due to 
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the small time step, we found it necessary to use 64-bit 
precision for these simulations. 

The simulation results are shown in Fig. 1131 where 
we see that there is only a small difference between 
the cooling performance in the adiabatic and full sim- 
ulations. Also shown is a simulation with a detuning 
A/27r = 2 GHz, evolved with substeps of At/2048, to 
increase the nonadiabatic effects. The parameters were 
scaled to compensate. Even with this relatively close 
detuning, we see that the cooling dynamics are well de- 
scribed in the adiabatic approximation. Other ensemble- 
averaged quantities corroborate this point. For the 4 
GHz detuning case, the atomic excited-state population 
(crV) is 8.4 X 10~^ (corresponding to a spontaneous emis- 
sion rate of 0.16); the adiabatic-approximation expres- 
sion g^|ap/ (A'^ + K^) predicts a similar value of 8.6x lO"'^ 
if we include a factor of 1 — (i?cff )/2Kiiax = 0.98 to ac- 
count for the ensemble-averaged spatial distribution of 
the atoms. For the 2 GHz case, the excited-state popu- 
lation is 3.4 X lO^'^ (spontaneous emission rate of 0.68), 
agreeing well with the adiabatic prediction of 3.5 x 10~^. 
In both cases, the cavity excitation (a^a) = 0.98, in good 
agreement with the adiabatic value of 1. 

VI. CONCLUSION 

Real-time quantum feedback control is worth investi- 
gating not only because of potential future applications, 
but because it will provide an undeniable litmus test of 
quantum measurement theory at the level of individual 
measurement records produced by measuring single sys- 
tems. In this work we have shown that, by using a simpli- 
fied estimation technique coupled with a relatively simple 
feedback algorithm, a single atom may be cooled close to 
its ground state in a realistic optical cavity. We have 
also explored the effectiveness of the feedback algorithm 
in various parameter regimes and we hope that this will 
help to guide future experiments. 
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APPENDIX A: SIMPLE FEEDBACK COOLING 
AND SQUEEZING 

Here we analyze the effect of the simple centroid-only 
feedback cooling algorithm described in Section IIII Bl 
This involves switching the optical potential to a high 
value when the atom is climbing a potential well, and 
switching it to a low value when it is falling. Recall that 
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the value of the "switched high" field is (1 + £i)E and 
the "switched low" field is (1 — £2)-^, so that the poten- 
tial height has the corresponding values {1 + ei)'^V and 
(1 — e2)'^V {E and V are the respective unmodulated 
values). 

Since the feedback algorithm involves more than sim- 
ply driving the system with an external force, to under- 
stand the action of the algorithm one must analyze its 
effect both on the means (the centroid) and the variances 
of the atomic position and momentum. The effect on the 
centroid is to damp the motion. In one cooling cycle, the 
energy of the centroid is multiplied by a "cooling factor," 

c•.=(i^)^ (Al) 

which is valid in the harmonic approximation. 

The effect on the variances is, however, to squeeze the 
atomic wave function. To understand why this is the 
case, consider what happens when the height of the po- 
tential is changed. This changes the frequency of the 
effective harmonic oscillator governing the atomic mo- 
tion, and as a consequence, a wave function that was not 
squeezed in the phase space of the old harmonic oscilla- 
tor, is squeezed in the space of the new one. In particu- 
lar, raising the height of the potential scales phase space 
so that momentum is squeezed, and conversely lowering 
the potential squeezes position. To determine the rate of 
squeezing, it is useful to consider what happens in one 
"cooling cycle." For example, we may define the cooling 
cycle to start when the atom crosses the bottom of the 
well, and end when the atom again crosses the bottom 
of the well traveling in the original direction (essentially 
one period of oscillation). Examining the effect of each 
of the changes in the potential during a cooling cycle, we 
find that the result of a single cycle is to multiply the 
momentum variance of the atomic wave function by a 
"magnification factor," 

«.^(i±|)', (A2, 

in the harmonic approximation. This factor is precisely 
the inverse of the cooling factor Ctb. 

This situation is analogous to an attempt to cool a 
classical ensemble of particles using a control algorithm 
designed for a single particle; although the algorithm will 
cool a single particle, it will almost certainly heat the 
other members of the ensemble. This will be an impor- 
tant consideration in designing an improved feedback al- 
gorithm. A simulation of this simple feedback algorithm 
is given in Figure |2 This shows that after a few cooling 
cycles in which the total energy is initially reduced, the 
resulting squeezing overtakes the cooling of the centroid, 
and the algorithm actually heats the atomic motion. 



APPENDIX B: HEATING MECHANISMS 

Spontaneous Emission. Spontaneous emission events 
kick the atom in both directions, and therefore cause mo- 
mentum diffusion. The effect of this is to increase the en- 
ergy, on average, linearly with time. This rate depends 
both upon the spontaneous-emission rate as well as the 
magnitude of the associated kicks. One point to note is 
that the magnitude of each kick is not merely determined 
by the momentum of the emitted photon, but also by 
the difference between the momentum of the excited and 
ground state wave functions. More precisely, we can fac- 
tor the atomic state vector just before a spontaneous- 
emission event into a product of internal and motional 
states, 

|V) = |V'c)|e) + |Vg)|g), (Bl) 

where \ipa) are states in the center-of-mass space of the 
atom, and "e" and "g" denote excited and ground atomic 
energy levels, respectively. Then in an unraveling of the 
master equation JQ) where the spontaneously emitted 
photons are detected but do not give position information 
about the atom, the one-dimensional atomic state just 
after a spontaneous-emission event is just the previous 
state with the internal state lowered and a momentum- 
recoil factor, 

|V) = |Vc)|g)e-'"^, (B2) 

where the random variable u is chosen from the projected 
angular distribution of the resonance fluorescence N{u) 
referred to in the SME iQ. Thus, in addition to the 
emission momentum recoil, the spontaneous emission ef- 
fectively converts the atomic state from the ground state 
wave function to the excited state wave function. In the 
regime where the adiabatic approximation is valid, we 
have the relation 

IV'e) = ^cos(fca:)|V'g), (B3) 

and thus we can interpret this extra process as being 
simply another momentum kick due to the absorption of 
a cavity photon (i.e., a superposition of a photon- recoil 
kick in either direction along the cavity axis). 

The actual rate of spontaneous emissions may be 
estimated by multiplying the spontaneous emission 
rate by the average excited state population, giving 
7g^|ap/(A^ + '*^) for an atom localized at a field antin- 
ode. Calculating this for the conditions that we use in 
the numerical simulations in the body of the text, we find 
that the rate is around 0.17. In the far-detuned regime 
that we study here, the spontaneous-emission heating is a 
minor effect compared to the other heating mechanisms. 
Measurement Backaction. The rate of heating from 
the back-action due to the continuous measurement pro- 
cess may be estimated by examining the effective master 
equation for the system in Eq. (|14|) . The energy of the 
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atomic motion in the effective potential, defined with re- 
spect to the minimum potential energy, is simply 



100 



l-cos\~kX) 



(B4) 



The rate of increase of this energy, averaged over all pos- 
sible trajectories, due to the measurement is determined 
by the term in the SME proportional to F, which gives 



= 87rr(cos2(fcX) sin^ikX)). 



(B5) 



Note that the T here is actually a time-averaged quantity 
when feedback is applied, which we assume to be well ap- 
proximated by the unmodulated value of F in Eq. (|17|) . at 
least to 0{e). The simplest limit for this expression is if 
we assume the particle to be randomly distributed along 
the potential, which is a reasonably good approximation 
to the initial heating rate: 



dt{E,B) 



= TrFfc^ 



(B6) 



A more appropriate limit for this heating rate after cool- 
ing has taken place is the low-temperature regime, where 
we can make the harmonic approximation so that the 
rate becomes 



dt{E, 



(B7) 



which is valid for small (Ecs)- 

Note that although we attribute this heating to 
measurement back-action, an equivalent, measurement- 
independent picture for this heating process is the cavity 
decay process. Because light escapes the cavity in dis- 
crete units at random times, there is a stochastic compo- 
nent to the cavity field intensity. Thus, there is a stochas- 
tic, heating force on the atoms. The rate obtained by this 
argument must be the same as the rate that we just ob- 
tained, since the heating rate obtained by averaging over 
all possible trajectories in a "quantum jump" unravel- 
ing of the master equation must be the same as in the 
"quantum diffusion" unraveling that we use here. 



APPENDIX C: EFFECTS OF DETECTION 
INEFFICIENCY AND FEEDBACK DELAYS 

1. Detection Efficiency 

Now we come to one of the two primary limitations in 
an experimental implementation of the present system, 
that of limited efficiency of the optical detectors. As this 
limits the information incorporated by the estimator, it 
is natural to expect that the best cooling is achieved for 
unit efficiency, and that the cooling simply gets worse for 
inefficient measurements. Simulated final energies using 
the Gaussian estimator are shown in Fig. 1141 which for 




FIG. 14: (Color online) Final energies, measured over the in- 
terval from f = 90 to 100, as the detection efficiency rj varies; 
other parameters match the canonical set. Cooling is accord- 
ing to the improved algorithm based on the Gaussian estima- 
tor; error bars reflect standard errors from averages over 128 
trajectories. 



the most part bear out this expectation. One remarkable 
feature of these results is that the cooling is extremely ro- 
bust to detector inefficiency. It is only down around 50% 
detection efficiency that a few trajectories begin to pull 
up the average energies, and below about 35% efficiency 
that there is a clear trend towards large steady-state tem- 
peratures. This robustness is very encouraging that an 
experimental demonstration of quantum feedback cool- 
ing will work well in spite of imperfect detectors. 



2. Speed and Feedback Delay Issues 

In the simulations we have discussed so far, the Gaus- 
sian estimator was updated every time interval of At^ = 
0.0005. For the canonical parameter set, this update in- 
terval corresponds to about 3 ns in physical time. This 
is obviously a very short time in which to evolve the es- 
timation equations H19|) as well as the curve fit outlined 
in Section Fill Bl Therefore, the first issue we will tackle 
is the speed with which one can iterate the control algo- 
rithm. 

Curve Fit. Recall that, given a set of estimates j/cst.n 
corresponding to times i„, we wish to implement the 
curve fit of the function oq + aix + a2X^ to the last q 
values of j/cst.i- The curve-fit coefficients are given by the 
solution of the normal equations. 
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where S„ -.^ and S'„ := ^2^=1 Vji^jT- How- 

ever, recomputing the sums in this system of equations at 
each iteration is not computationally efficient. To make 
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the on-the-fly calculation of this curve fit feasible, we can 
instead implement the coupled recurrence relations 
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and then compute the fitted slope coefficient at time t,- 
according to 
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(C3) 

For notational convenience, we have defined the times 
corresponding to the last q measurements to be i„ = 
1 — n, so that the fitted slope at the current time is simply 
ai. Thus, we can simply trigger the feedback cooling 
algorithm on the sign of oi, switching the potential high 
or low for a positive or negative sign of oi , respectively. 
The speed of this algorithm is largely independent of the 
number of samples used in the curve fit, assuming that 
the table of the last q values can be accessed quickly (i.e., 
it can fit into fast memory). 

Later we will examine the effects of computation and 
signal propagation delays on the cooling performance. 
Such delays are easily accounted for in the present al- 
gorithm by using the curve fit to extrapolate forward in 
time. To do this, we require the quadratic fit coefficient. 
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so that we trigger on the value of the current-time slope, 
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(C5) 



where d is the delay in multiples of Ate from the time of 
the most recent data used in the fit to the current time. 
Due to the predictive nature of this fitting scheme, we 
expect that the cooling should be robust to small delays. 

Gaussian Estimator. The other issue to be resolved 
is to maximize the efficiency of updating the Gaussian 
estimator. One computationally expensive aspect of the 
form of the equations H19|l is the need to evaluate tran- 
scendental functions. A substantial speed improvement 
results from introducing the variables Xi := sm(2k{X)c), 
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cos(2fc(X)e), (P)c/fc, Xi := 2PV^, 



exp(-2A:2y^), xq := Vji/P, and := exp(-2PF^), we 



can rewrite the estimator equations in the form 
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dxe = 

dxy = 
with 



(47rPx2a;3 - ArjTxiX^x^jdt 
+ \/8r/TxiX2X4X5dWe 

— {■iTTpxiXs + 'ir]rxfx2X^x^)dt 
—^/SrJTxiXiXzdWc 

— VmaxXlX^dt + ^/SrJTxiX^XTdWo 

{SnPxj - 4:i]rx^xlx^)dt 
+^/8r]rx2xix5dWc 

[Snk'^x^xj + 4r]T{xiX5 — X2xl)x'lx^]dt 

— \/?>rjTx2 x^x^dWc 
("4V;„ax2:2a:5a;7 + r[l + ^^'(l - 2a;|)] 
-8r]Txlxix^)dt 

-y2^[l - 2{x4 + 2x^)xz]x2dWe 
(27rfc^X6 — Vmaxa;2a;4a;5 — 'i-qVxiXix1x7)dt 
+ -y8r/rx2 a;4 X5 xy d We , 



(C6) 



dW, = 



'2r]r{l + X2Xz)dt, 



(C7) 



thus obviating any need for transcendental- function eval- 
uation. However, it turns out that the evolution of these 
equations via the Euler method is much less stable than 
the evolution of the original equations. This is due to 
the fact that there are pairs of variables [{xi^X2) and 
[x^^^xc,)] that evolve separately but must satisfy consis- 
tency conditions [x\ + X2 — ^ and X5 — exp(— X4)]; also 
a large stochastic fluctuation can push these values into 
unphysical ranges (e.g., xi > 1). A better strategy is to 
emulate the results of an Euler solution to the original 
estimator equations, and hence update xi, X2, and 0:5 
using the increments 



X2 sm{2kA{X)c) 

xi sin(2A:A(X)o) (C8) 



Axi = j:i[cos(2fcA(X)c) - 1] 
Ax2 = a;2[cos(2fcA(X>c) - 1] 
Ax5 — X5 exp(— Aa;4), 

where 

2kA{X)c ^ ink'^xsAto + ^/8^xiX4X5AWc, (C9) 

AWait) := //'^^*'"dWc, and Ax^ is the Euler-update 
increment from Eqs. HC6|I . These equations can then be 
expanded in powers of 2kA{X)c to avoid computing the 
transcendental functions. A final trick to help stabilize 
the equations is to multiply the xi and X2 variables at 
each time step by the normalization quantity [1 — {xi + 
)]/2, which explicitly maintains the identity x\-\-x\ = 
1 without requiring the evaluation of a square root. 
Execution Times. To evaluate how quickly the control 
algorithm can be iterated, we implemented the control 
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FIG. 15: (Color online) Final energies, measured over the 
interval from t = 90 to 100, as a function of the feedback 
delay rj. (Time is measured in units of the 5.8 fis harmonic- 
oscillator period.) Cooling is according to the improved al- 
gorithm based on the Gaussian estimator, with extrapolation 
implemented in the curve fit; error bars reflect standard errors 
from averages over 128 trajectories. Other parameters are as 
in the canonical set. 



algorithm on a modern, general-purpose, serial micro- 
processor (1.25 GHz Alpha EV68). A single iteration 
of the estimator according to Eqs. H19() and a curve fit, 
including all steps necessary to compute the control out- 
put, can be completed in about 260 ns. Switching to 
the faster method of Eqs. (|C6I) with the more stable up- 
dates (jCSII can be completed in 140 ns with a truncation 
at fourth order in 2kA{X)c, and 120 ns with a second- 
order truncation (including a normalization step). The 
cooling performance is essentially the same with all these 
algorithms, in that they all reach the final temperature; 
the exception is an algorithm that uses Eqs. (|C6p with- 
out the modifications of Eqs. ljC8|l . which has substan- 
tially worse cooling performance. The tradeoff in the 
truncated-expansion algorithms is that lower-order trun- 
cations are faster but also more unstable and thus require 
more resets of the estimator. 

Based on this performance, and assuming Moore's 
law, general-purpose hardware will not be able to iter- 
ate this algorithm in 3 ns for about a decade. However, 
even presently available, specialized technologies such as 
multicore processors and field-programmable gate arrays 
(FPGAs) can likely be used to significantly decrease the 
execution time. For example, because the curve fit can 
be moved onto a separate logical processor, the estimator 
can be updated more quickly (84 ns for the fourth-order 
version, 70 ns for the second-order version); the extra 60 
ns required to generate the feedback value via the curve 
fit then represents an additional delay, but does not limit 
the rate at which the estimator can be iterated. Further 
parallelization should bring these times down substan- 
tially. 

Impact on Cooling Performance. Given these lim- 
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FIG. 16: (Color online) Final energies, measured over the in- 
terval from t = 90 to 100, as a function of the Gaussian esti- 
mator time step Ate- (Time is measured in units of the 5.8 fis 
harmonic-oscillator period.) Shown are results obtained by re- 
ducing the number of samples to keep the time interval for the 
quadratic curve constant (circles) as well as results obtained 
via the multiple-estimator staggering method to maintain a 
300-point fit (triangles). Error bars reflect standard errors 
from averages over 128 trajectories. Other parameters are as 
in the canonical set. 
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FIG. 17: (Color online) Ensemble energy evolution illustrat- 
ing the combined effects of detection efficiency -q, feedback 
delay Td, and estimator time step Ate- (Time is measured in 
units of the 5.8 fis harmonic-oscillator period.) These param- 
eter values for the curves (a)-(d) are given in the text; others 
are as in the canonical set. Extrapolation is implemented 
in the curve fit, and multiple-estimator staggering is used to 
maintain a 300-point flt. Each curve represents an average 
over 128 trajectories. 



itations in computing power it is necessary to evaluate 
how quickly the algorithm needs to be iterated in order 
to achieve good cooling performance. The cooling per- 
formance in the presence of feedback delay = dAt^ is 
shown in Fig.^| The cooling is very robust to the delay, 
since we can compensate by extrapolating the curve fit. 
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The cooling behavior is essentially unaffected out to a de- 
lay of about 460 ns, which is well above the 60 ns figure 
that we mentioned for generating the feedback value, and 
thus additional signal propagation and processing delays 
should not be too problematic. 

Limitations on the time Ate required to iterate the 
Gaussian estimator can also seriously impact the cool- 
ing performance. Not only does a large Ate step re- 
duce the quality of the estimated solution, but in fit- 
ting the quadratic curve over a finite interval, there are 
fewer points involved in the fit, reducing the effectiveness 
of the fit's noise reduction. These effects on the cool- 
ing behavior are illustrated in Fig. ^1 which shows that 
the late-time ensemble energies grow approximately lin- 
early with Ate- Here, the 70 ns figure mentioned above 
for Aio yields a final energy of around 25, compared to 
around 9 for the smallest value of At^ = 0.0005 (3 ns) 
shown, a substantial decrease in the cooling effectiveness. 
One possible way to address this problem is to use multi- 
ple Gaussian estimators that are staggered in time. For 
example, defining Atj^in = 0.0005, and given that the 
computational hardware is limited to A^q = NAt^iin 
for some integer N, we can arrange to have N inde- 
pendent implementations of the Gaussian estimator ex- 
ecuting in parallel. They all receive the same measure- 
ment record as input, but they are staggered in time 
such that one produces output at each small time step 
Ainiin — 0.0005. In this way, each of the individual so- 
lutions is of lower quality than a fine-time-step solution, 
but the curve fit can be fit to all of them (restoring a full 



300 point curve fit), thus helping to stabilize the cooling 
algorithm. The implementation of this method is also 
displayed in Fig. ^| which shows that this method re- 
stores much of the quality of the cooling that was lost in 
using the simpler algorithm. The rightmost point in this 
figure corresponds to operating 55 parallel estimators. 

So far, we have studied the performance impact on 
cooling of each of the three equipment-limitation effects 
(detection efficiency, feedback delay, and estimator time 
step) separately. Some representative examples of their 
combined effect are shown in Fig. 1171 The four curves, 
in order of increasing departure from the ideal case, cor- 
respond to the following parameters: (a) r] — 1, — 0, 
Ato — Aimin = 0.0005 (i.e., 3 ns in physical time); (b) 
ri = 0.8, Td Ato = lOAimin = 0.005 (i.e., 29 ns); 
(c) 7] = 0.7, Td = Ato = 25Aimin = 0.0125 (i.e., 72 
ns); (d) 77 = 0.5, Td = At,, = 40Ai„iin = 0.02 (i.e., 116 
ns); (e) ry = 0.3, Td = Ata = 50Atmin = 0.025 (i.e., 
145 ns). The sums of the populations in the lowest two 
bands for these cases, computed between t = 90 and 
100, are (94 ± 5)%, (84 ± 5)%, (45 ± 4)%, (47 ± 4)%, 
and (13 ± 2)%, respectively. The cooling performance, 
even for the case that should be within reach of current 
technology (d), is not unreasonable in comparison to the 
ideal case (a), although cooling quickly becomes worse 
as these parameter values are relaxed (e). Even with 
these effects, this method is certainly suitable for long- 
time storage of atoms in the optical potential and can 
still reach the ground state with serviceable efficiency. 
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